#!/bin/bash
# (1) Enumerate different isoforms for a "gene" based on Tophat's output. 
# (2) Call PerM to map reads on isoforms. Generate *.mapping and *.isoforms.
# (3) Call TranSEQ.sh to generate the input for TranSEQ.
# In TranSEQ.sh, the abundance level of the enumerated isoforms are estimated.
# The required scripts should be in ~/SCRIPT/ and required shells should be in ~/SHELL/
# PerM and TranSEQ are set in $PATH by .bashrc
PANF=/home/cmbpanfs-01/yanghoch #can be overwritten by bashrc
source ~/.bashrc
echo $0
if [ $# -ne 2 ]; then
	echo "The 1st arg is the refName."
	echo "refName.fa, refName.sea and ref.bed are used to enumerate isoforms"
        echo "The 2nd arg is the reads path, used to evaluate the abundance level"
	echo "Two arguments are needed."
	echo "The out is"
	echo "(i) Expression levels and LR in \$1.out."
	echo "(ii) Selected isoforms in $1_CDSs_selected.fa"
	exit
fi
READS=$2
echo ""
echo "(i) Enumerate all possible isoforms with long enough CDS "
if [ ! -f $1.fa ]; then
	if test -s ../$1.fa; then
		ln -s ../$1.fa $1.fa
	else
		echo "$1.fa not found."
		exit
	fi
fi
$PANF/EnumI $1
RefSEQ=$(basename $1)'_I'
cp $1'_I.fa' $RefSEQ'.fa'
cp $1'_CDSs.fa'  $(basename $1)'_CDSs_selected.fa'
RefPrefix='I'
#RefSEQ=$(basename $1)_CDSs
#mv $1_CDSs.fa $RefSEQ.fa
#Temp use all CDSs as slectedted one
#RefPrefix=CDS
if test ! -s "$RefSEQ.fa"; then 
	echo "No isoform is generated or passed the filtering"
	exit
fi
echo "Get the isoform names and lengthes " 
python ~/SCRIPT/one_AceView_gene_one_isoforms_file.py\
 $RefSEQ.fa $RefPrefix $1.isoforms
echo "Get the first isoform sequence"
#get the first isoform for plot 
bFlag=T
~/SHELL/getSeqsFromFasta.sh $1'_I.fa' I $bFlag I0.fa
if [ -f $1.isoforms ]; then 
	cp $1.isoforms $(basename $1).isoforms #save a copy 
else
	echo "$1.isoform is not created!!"
	exit
fi

# A naive cut to limit the size of isoforms
#head -n 10000 RefSEQ.fa > tmp.fa
#mv tmp.fa RefSEQ.fa
#echo "(ii) Map reads to the enumerate isoforms"
echo ""
echo "(ii) Map reads to the enumerate coding seqeunces "
$PANF/PerM $RefSEQ'.fa' $READS 5 -B -m > PerM_$(basename $1).log
\rm PerM_$(basename $1).log
\rm *miss.fasta 2> tmp.log
\rm tmp.log
# concatenate all mapping result to $1.mapping
baseReadFileName=$(basename $READS ".txt")
~/SHELL/catFiles.sh ".mapping" $baseReadFileName $1.mapping

if test ! -s "$1.mapping" ; then
	echo "$1.mapping not found" 
	exit
fi
echo ""
echo "(iii) Call TranSEQ.sh to select the expressed isoform"
TranSEQ_Out=$(basename $1).out
~/SHELL/TranSEQ.sh $1 $TranSEQ_Out 
python ~/SCRIPT/select_isoforms_from_TranSEQ_out.py\
 $RefSEQ.fa $TranSEQ_Out $RefSEQ'_selected.fa' 

